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ABSTRACT 

One of the tenets of the radio pulsar observational picture is that the integrated 
pulse profiles are constant with time. This assumption underpins much of the fantas- 
tic science made possible via pulsar timing. Over the past few years, however, this 
assumption has come under question with a number of pulsars showing pulse shape 
changes on a range of timescales. Here, we show the dramatic appearance of a bright 
component in the pulse profile of PSR J0738— 4042 (B0736— 40). The component arises 
on the leading edge of the profile. It was not present in 2004 but strongly present in 
2006 and all observations thereafter. A subsequent search through the literature shows 
the additional component varies in flux density over timescales of decades. We show 
that the polarization properties of the transient component are consistent with the 
picture of competing orthogonal polarization modes. Faced with the general problem 
of identifying and characterising average profile changes, we outline and apply a sta- 
tistical technique based on a Hidden Markov Model. The value of this technique is 
established through simulations, and is shown to work successfully in the case of low 
signal-to-noise profiles. 
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1 INTRODUCTION 

The radio emission from pulsars is characterised by a range 
of dynamic phenomena that take place on various timescales. 
Microstructure is observed at the shortest (/xs) timescales, 
stochastic or organised changes in the pulse shape occur on 
timescales of the rotational period P, while other phenom- 
ena such as nulling (e.g. Lorimer & Kramer 2005), where 
radio emission totally switches off, may last significantly 
longer than one rotation. An examination of the average 
pulse properties of a given pulsar however, demonstrates 
that the mean shape of a sufficiently large number of in- 
dividual pulses remains remarkably constant over long pe- 
riods of time. It is this quality of radio pulsars that makes 
them extremely useful tools: knowing the exact shape of the 
pulse profile increases the precision of the measurement of 
the pulse time of arrival, which is the basic quantity in pulsar 
timing experiments. 

Although it is rare that the integrated profiles of pul- 
sars change with time, it is not unprecedented. Timescales of 
changes range from hours through to months and years. On 
short timescales are the rotating radio transients (RRATs) 



(McLaughlin et al. 2006), which are thought to be neutron 
stars that only emit individual bursts of emission at irregu- 
lar and infrequent intervals. A group of pulsars show mode 
changing, where two distinct and different pulse profiles are 
observed over timescales of hours (e.g. Gil et al. 1994). The 
intermittent pulsar PSR B1931+24 (Kramer et al. 2006) is 
present for 5 to 10 days before its emission ceases for some 
30 days. In this case, the derivative of the rotational period, 
P, changes between the on and off phases and is therefore 
correlated with changes in the radio pulse profile. Similar 
results associated with more subtle profile changes are re- 
ported in Lyne et al. (2010), and another intermittent pulsar 
(PSR J1832+0029) is discussed in Kramer et al. (2008). In 
PSR J1119-6127 pulse changes were seen which appeared to 
be associated with glitch activity (Weltevrede et al. 2011). 
There is strong evidence that all these effects (nulling, mode 
changing and intermittency) are magnetospheric in origin. 

A further long-term effect, due to geodetic precession 
of a pulsar in a binary orbit, can also cause pulse shape 
changes. These have been observed in e.g. PSR B1913+16 
(Kramer 1998, Weisberg & Taylor 2002) and Jl 141-6545 
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Figure 2. A greyscale representation of the intensity of 71 pro- 
files, observed at irregular intervals between 2003 and 2011. The 
peak intensity of each profile is normalised to unity; the pro- 
files are aligned by cross-correlation with a top-hat function. The 
change occurs at profile 23. 



(Manchester et al. 2010). Periodic changes in the average 
profile of PSR B1828— 11 were interpreted as free precession 
by Stairs et al. (2000), however free precession of solitary 
neutron stars is considered unlikely (Sedrakian et al. 1999) 
on theoretical grounds. Changes in the pulse profile due to 
precession relate to changes in the viewing geometry rather 
than magnetospheric effects. 

The average pulse profiles of pulsars are partially lin- 
early polarized, to a lesser or higher degree. Highly energetic 
pulsars feature high degrees of linear polarization (e.g. Wel- 
tevrede & Johnston 2008). There is good evidence to suggest 
that the observed emission results from the superposition of 
two orthogonally polarized modes (OPM), which arise and 
propagate inside the pulsar magnetosphere (e.g. McKinnon 
& Stinebring 2000, Karastergiou et al. 2001). Comparable 
intensities of the modes has been favoured observationally 



as the cause for reduced linear polarization in pulsars at 
typical (~1 GHz) observing frequencies. A further observa- 
tional consequence is that the changes in the structure of 
the total power average profile with observing frequency are 
often coupled with particular changes in the degree of linear 
polarization, reflecting the spectral behaviour of the orthog- 
onal polarization modes (as discussed in Karastergiou et al. 
2005, Smits et al. 2006);as the OPMs become more equal 
in strength, the total power increases and the polarization 
decreases. 

Pulsar timing models need to incorporate all known 
physical phenomena (intrinsic to the pulsar or not) that af- 
fect the measured times-of-arrival, in order to achieve a floor 
of sensitivity that would enable the discovery of extremely 
weak components to the model, such as gravitational waves 
(e.g. Hobbs et al. 2009). Pulsar timing uses template match- 
ing and relies on the average profile not varying with time. 
Any variability in the profile adversely affect the timing 
model. In the following, we present data from a significant 
change in the average pulse profile of PSR J0738— 4042. We 
show how polarization data reveal details about the change, 
and discuss possible interpretations. We present a robust 
statistical technique to characterise the change, and explore 
its potential physical origins. 



2 THE TOTAL INTENSITY PROFILE OF PSR 
J0738-4042 

PSR J0738— 4042 was one of the first radio pulsars discov- 
ered (Large et al. 1968). It has a high dispersion measure of 
160.8 cm~ 3 pc, and its spin period of P = 375 ms and a pe- 
riod derivative of P = 1.61 x 10~ 15 places it within the bulk 
of normal pulsars on the P-P diagram. Its relatively high 
flux density (80 mjy at 1.4 GHz) has made it a consistent 
observing target over the 40 years since its discovery. The 
most recent polarimetric profiles of the pulsar over a wide 
range of frequencies have been published in Karastergiou & 
Johnston (2006), Johnston et al. (2006) and Johnston et al. 
(2007). Average profiles in three separate observing bands 
taken in 2004 are shown in Fig. [I] (thin line) . The profile at 
1.4 GHz shows one bright component with a shoulder com- 
ponent on its leading edge, on what appears to be a broad 
pedestal of emission on the leading and trailing edge. Sub- 
sequently, the pulsar was observed as part of a program to 
measure accurate rotation measures (Noutsos et al. 2009). 
The 1.4 GHz profile, taken in 2006, is shown as the thick line 
in Fig. [I] An additional component, ~ 15° earlier than the 
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Table 1. 40 years of average profiles from PSR J0738-4042. 



Date 


Frequency 


Component at -15° 


Reference 




<1970 
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<1975 


1400 MHz 


Strong and discrete 


Backer (1976) 




<1977 


631 MHz 


Shoulder to main pulse 


McCulloch et al.(1978) 




<1977 


1612 MHz 


Strong and discrete 


Manchester et al. (1980) 




1979 


950 MHz 


Strong and discrete 


van Ommen et al. (1997) 




1990 


950 MHz 


Weak shoulder to main pulse 


van Ommen et al. (1997) 




1991 


800 MHz 


Absent 


van Ommen et al. (1997) 




1996 


1375 MHz 


Absent 


unpublished 




1997 


1375 MHz 


Absent 


unpublished 




2004 


1375 MHz 


Absent 


Karastergiou & Johnston 


(2006) 


2004 


3100 MHz 


Absent 


Karastergiou & Johnston 


(2006) 


2005 


8400 MHz 


Absent 


Johnston et al. (2006) 




2005 


3100 MHz 


Absent 


Johnston et al. (2007) 




2006 


1369 MHz 


Strong and discrete 


Noutsos et al. (2009) 





main peak, can be seen on the leading edge of the 2006 pro- 
file, which is almost entirely absent in 2004. This component 
is present at all three observing frequencies. 

Armed with this result, we looked through the literature 
for other published profiles of this pulsar. Table [l] includes 
a summary of available data, where the date and observing 
frequency are given with a note relating to the presence of 
the leading component. Two facts are immediately evident. 
First, there is a period between 1991 and 2005 where the 
component is totally absent. Secondly, the presence or ab- 
sence of the leading component is a broadband phenomenon; 
there are no contradictory observations at a particular fre- 
quency. This is clearly seen in the representative, high S/N 
profiles in Fig. [I] all showing additional leading edge emis- 
sion in 2006. 

Figure [2] shows data from 71 pulse profiles taken with 
the 21 cm receiver of the Parkes telescope. The profiles, 
which are not collected at regular intervals, span a range 
of dates from February 2003 to January 2011. The index 
of each profile as shown on the y-axis is used for the sta- 
tistical analysis that follows. The duration of each observa- 
tion varies from 120 to 1199 s. All profiles in Fig. [2] have 
been normalised to the peak intensity and aligned by cross 
correlation with a top hat function. Although absolute flux 
density calibration is not available for all observations, an 
analysis of the S/N as a function of observing time reveals 
that the peak flux density remains broadly similar, within 
a ~ 20% margin. The crucial change in the shape of the 
leading edge can be clearly detected after profile #23 with 
the "new" component remaining on all subsequent profiles. 
Towards the end of the sequence the new component has 
grown to its brightest amplitude. 
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3 CHANGES IN THE POLARIZATION 
PROFILE 

As mentioned in the introduction, polarization is a useful 
diagnostic in the interpretation of pulsar radio emissions. 
Figure [3] shows three average profiles of PSR J0738-4042, 
from 2004, 2006 and 2010. It is immediately obvious that 
the additional component which appears between 2004 and 
2006 is also associated with a different polarization state. 



Figure 3. Three average profiles of PSR J0738-4042 in full po- 
larization, from 2004, 2006 and 2010 (top to bottom). The total 
intensity (black), linear polarization (red) and circular polariza- 
tion (blue) are shown just below the PA curve (dotted line). The 
extra component in the 2006 and 2010 profiles is clearly respon- 
sible for the observed additional orthogonal jumps in the leading 
edge of that profile. These jumps coincide in phase with local 
minima in the linear polarization. 
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Figure 4. Simulated polarization profiles resembling the data of 
2006 and 2010 from Fig. [3] These profiles have been generated by 
adding an orthogonally polarized, Gaussian component of fixed 
width and variable amplitude to the polarization profile without 
the transient component (from 2004). Most features of Fig. [3] 
around pulse phase —15° are well reproduced. 



Between pulse phase —20° and —10°, there are 3 orthog- 
onal polarization jumps in the 2006 and 2010 profiles, as 
opposed to a single jump in the 2004 data. It is also evi- 
dent that the degree of linear polarization in this region of 
the profile is related to the total power, and that there are 
local minima directly related to the orthogonal PA jumps. 
A comparison between the 2004 and 2006 data shows that 
as the total power increases around pulse phase —16°, the 
linear polarization drops. 

We attempted to reproduce the observed changes in po- 
larization using a simple model. The starting point of the 
model is the polarization profile obtained from observations 
in 2004 shown at the top of Fig. [3] To this, we added a Gaus- 
sian component centred at phase —15.47°, which is 100% po- 
larized. We computed the total intensity I c of the simulated 
component as: 



I c (0) = Aexp [-(</» - 15.47) 2 /(2<7 2 )] 



(1) 



where A the amplitude and a the width of the Gaussian, 
and <f> the pulse phase. For each pulse phase bin, we set the 
polarization of the simulated component to be orthogonal to 
the polarization of the 2004 profile, by computing its Stokes 
parameters (I C ,Q C ,U C ,V C ) relative to the Stokes parameters 
of the 2004 profile (Io,Qo,U ,V ), taking into account that: 

J^ = ^,ic = v^!+W+v! (2) 

*^c 

Therefore, Q c , U c and V c can be computed for every pulse 
phase bin as: 

Qolc 



Q c = 



(3) 



U c 



Uolc 



VQ 2 o + ul + vg 



(4) 



(5) 



where the — signs ensure that the polarization vectors are 
antiparallel on the Poincare sphere. Finally, we add the 
Stokes parameters of the new component to the 2004 data 
to produce simulated polarization profiles. 

Figure [4] shows the best results of our model, which 
can be compared directly to Fig. [3] We find that we can 
reproduce the data very closely by setting a — 1.43° and 
varying A between ~0.25 for the 2006 data and ~0.5 for the 
2010 data. The fact that, at the relevant pulse phase region, 
the 2004 profile is almost entirely polarized, and the 2006 
and 2010 profiles can be simulated by adding a 100% and 
orthogonally polarized component leads to the conclusion 
that the polarization profiles with the additional component 
are direct observations of the superposition of two, 100% 
polarized, orthogonal modes of emission. 



4 STATISTICAL DESCRIPTION OF THE 
PROFILE EVOLUTION 

4.1 Hidden Markov models 

Visual inspection of the profiles of Fig. [2] suggests that the 
profile evolution of PSR J0738— 4042 can be described by 
a single change of state, an assertion that is supported by 
the polarization modelling described in the previous sec- 
tion. We utilize a fully probabilistic realization of a hid- 
den Markov model (HMM) to automatically identify pu- 
tative state changes in the data. Hidden Markov models 
(Rabincr 1989) have been widely used for inferring latent 
changes in sequential data. Consider a sequence of observa- 
tions, Y = {y t }f =1 , y t G R d ,Vi. The distribution of an ob- 
servation, yt is determined by a corresponding hidden state, 
St £ {1, . . . , J} and a state dependent observation probabil- 
ity. 

The hidden state at time t = 1 is determined by a prior 
state vector, tv = [tti, . . . ,7rj] T , where Hj = p(si = j). The 
Markov property implies that a hidden state st-i, s t depends 
only on st-i and not on those at time t — 2 and before: 



p(s t \s t -i, ■ ■ ■ , si) = p(s t \s t -i). 



(6) 



State transitions are jointly determined by a transition ma- 
trix, A = [aij], where mj = p(st ~ j\st-i ~ i) and obser- 
vation likelihoods. An observation yt depends only on the 
corresponding hidden state, St and we utilize here a state- 
dependent Gaussian density, such that the predictive distri- 
bution over observations, conditioned on the hidden state, 
is given by: 



P(yt\st =j)=N (yt; fip Sj) 



(7) 



where /1 ■ and Ej are a the mean vector and a covariance 
matrix for the j th state. 

Inference in a HMM proceeds by evaluation of the pos- 
terior distribution over the parameters, {tt, A, fi, S} as well 
as the posterior distributions of the hidden state variables, 
p(s t \Y). A two-stage maximum-likelihood algorithm is often 
adopted, such as the Baum- Welch algorithm (Baum et al 



19701, a special case of the expectation-maximisation (EM) 
algorithm ([Dempst er et al.|1977[ ). The most probable state- 
sequence can be found using the Viterbi algorithm ( Rabiner 
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Figure 5. Contours of the entropy of the HMM, superposed on 
the profiles from Fig. [2] The entropy is a measure of the prob- 
ability of a state change, and reaches the value 1 in profile 23 
and pulse phase -16. The contours are at entropy value steps of 
0.1, from to 1. Apart from profile 23, the probability of a state 
change is increased in profile 68. 



1989). The maximum-likelihood approach, however, has ma- 
jor limitations; most notably overfitting and inference of un- 
derlying model complexity, such as determining the most 
probable number of states. In order to address these limita- 
tions, we employ a fully Bayesian approach, which exploits 
the tractable bounds of variational Bayes approximations 
( Jordan et al.| |T999). The Bayesian methodology allows for 
uncertainty to be handled at all levels of inference, mak- 
ing the approach ideal for analysis of smaller data samples. 
Furthermore, the number of underlying states is inferred 
automatically such that sequences without significant state 
changes are modeled by a single state process. This deviates 
significantly from maximum-likelihood methods in which the 
data is forced into a set number of states. 

As posterior state probabilities are inferred for each da- 
tum, we can investigate with ease not only the existence 
of state changes but also track the entropy associated with 
state determination at each point. The information entropy 
associated with the posterior over the states, conditioned on 
observation y t is hence given as: 



Ht = ~^P(st = j\y t ) logp(s t =j\yt)- 



(8) 



The advantage of utilizing entropy lies in the fact that in- 
creases in entropy will occur in locations in which a full 
state-change is not supported by the data. If the logarithm 
is to base two, then the entropy is returned in bits, with 1 
bit of entropy indicating the existence of a fully supported 
state change. 

We use the Bayesian HMM to infer the posterior dis- 
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Figure 6. Results of the sensitivity analysis. The curve shows the 
fraction of 1000 simulated profiles where the HMM successfully 
identifies a state change, versus the S/N of the new component. 
The HMM performs perfectly for S/N of 1 and above. 



tribution over the state sequence for each bin of the pulse 
profile, using the total power amplitudes of each bin as our 
observables. Figure [5] shows the entropy (in bits) of the state 
posterior applied to all bins of all profiles, in contours su- 
perposed on the raw data. The contours, which range from 
to 1 in steps of 0.2, indicate very high probability of a 
state change in profile 23, but also suggest there may be 
further change in the most recent observations (profile 68). 
Artificially changing the order of the profiles moves the iden- 
tification point of the state change accordingly. 



4.2 Sensitivity analysis 

Although the HMM clearly and robustly identifies a state 
change in the average pulse profile, it is not surprising that 
this analysis performs well, given the large magnitude of 
the event (it can be picked out clearly by eye). We have per- 
formed a sensitivity analysis to test the performance of the 
HMM on noisier data. This involves adding increasing levels 
of white noise to the profiles and comparing the state-change 
identification with the original analysis. We have performed 
this analysis by means of a Monte Carlo simulation, pro- 
ducing 1000 versions of each profile with a given S/N and 
changing the noise from 4 to 100 times the original in steps 
of 4 (i.e. 25 x 10 3 simulated profiles for each observed pro- 
file). For each set of profiles at a given S/N, we run the 
HMM and count the fraction of sets in which a state change 
is detected, as well as the profile number that the transition 
occurs. 

Figure [6] shows the results of the sensitivity analysis, 
as the fraction of the 1000 sets where a state change is de- 
tected, versus the S/N of the peak of the new component in 
the simulated data. Addition of the noise levels mentioned 
above artificially decreases the average S/N of the peak of 
the new component to values between 0.4 and 10. The HMM 
performs extremely well from S/N per pulse phase bin of 1 
and above. 
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5 DISCUSSION AND CONCLUSIONS 

We have shown that PSR J0738-4042 has undergone dra- 
matic changes in its average pulse profile in the period be- 
tween 2004 and 2006. The changes affect the total amplitude 
and polarization of a single component on the leading edge 
of the profile and are broadband at least over the frequency 
range between 600 and 3100 MHz. Examination of the lit- 
erature shows that the transient component was present be- 
tween 1970 and 1990, then absent until 2006 since when it 
has again been present up to the current epoch (January 
2011). The 'intermittency' of this particular component can 
therefore be measured in tens of years compared to e.g. the 
tens of days in PSR B1931+24. This serves as a unique ex- 
ample of magnetospheric changes on very long time scales. 
This is further proof that there remains a lot to be under- 
stood on the dynamic nature of pulsar magnetospheres on 
all timescales. 

The increase in intensity of the transient component is 
accompanied by an orthogonal transition in the position an- 
gle of the linear polarization, and by a decrease in the frac- 
tional linear polarization of that part of the profile. When 
the extra component is absent the leading part of the profile 
is almost completely polarized. Orthogonal jumps in the po- 
larization angle are common in pulsars with medium or low 
levels of linear polarization, which implies that the degree 
of polarization is affected by the presence of the two modes. 
In the past, two models have been put forward to account 
for this behaviour: the observed pulsar radiation either oc- 
curs in two partially polarized orthogonal modes of emis- 
sion which are emitted disjointly (e.g. Cordes et al. 1978), 
or two entirely polarized orthogonal modes, the superposi- 
tion of which sets the total polarization (eg McKinnon & 
Stinebring 2000). The data presented here strongly favour 
the latter, as this is the first example where we see two dis- 
tinct states: a single mode and high polarization when the 
additional component is absent, and low polarization asso- 
ciated with higher total intensity when it appears. We show 
the validity of this interpretation of the post-2006 data with 
a simple simulation. 

Orthogonal polarization modes are thought to be re- 
lated to propagation effects within the pulsar magnetosphere 
(e.g. Melrose 2000), which then strongly suggests a magneto- 
spheric rather than geometric origin for the observed profile 
changes. We consider a geometrical effect such as free preces- 
sion to be extremely unlikely. Magnetospheric effects have 
been shown to be responsible for the profile changes seen 
in PSR B1931+24, and the correlations between the period 
derivative and the profile shape changes of a small number 
of pulsars in Lyne et al. (2010) also point to magnetospheric 
effects. In PSR J1119-6127, profile changes seem to occur 
immediately following a glitch (Weltevrede et al. 2010). Al- 
though there are insufficient existing data to trace the be- 
haviour of the timing of the pulsar over the long term, there 
is no evidence for glitch activity between 2004 and 2006. 
Also, unlike PSR J1119-6127, a long series of single pulses 
from PSR J0738-4042 taken after 2006 shows the transient 
component to be persistent in the single pulses with no more 
than typical variability. 

We have presented a robust statistical technique based 
on a Hidden Markov Model that estimates the likelihood of 
state changes in the total power data, and showed that this 



technique has great potential for substantially noisier data 
in which similar systematic changes occur. We plan to apply 
the HMM technique to a large number of pulsars for which 
regular observations have been conducted, to look for and 
characterise statistically significant profile variations. Most 
importantly, this technique provides a means of correlating 
changes in the pulse profile with other physical parameters, 
such as the period and period derivative, which could prove 
extremely useful in partially accounting for non-Gaussian 
timing residuals in pulsar timing models. 
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